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Abstract: In this paper we introduce and analyze a class of diffusion type equations related to certain non- 
Markovian stochastic processes. We start from the forward drift equation which is made non-local in time by the 
introduction of a suitable chosen memory kernel K(t). The resulting non-Markovian equation can be interpreted 
in a natural way as the evolution equation of the marginal density function of a random time process l(t). We 
then consider the subordinated process Y(t) = X(l(t)) where X(t) is a Markovian diffusion. The corresponding 
time evolution of the marginal density function of Y(t) is governed by a non-Markovian Fokker-Planck equation 
which involves the memory kernel K(t). We develop several applications and derive the exact solutions. We 
consider different stochastic models for the given equations providing path simulations. 

1 Introduction 

In this introduction, we describe and motivate the themes developed in the paper. Historical notes will be 
presented in Section [2j 

Brownian motion B(t), t > 0, is a stochastic process with many properties. It is at the same time Gaussian 
and Markovian, has stationary increments and is self-similar. A process X(t), t > 0, is said to be self-similar 
with self-similarity exponent H if, for all a > 0, the processes X(at), t>0, and a H X(t), t>0, have the same 
finite-dimensional distributions. Brownian motion is self-similar with exponent H = 1/2. In contrast, frac- 
tional Brownian motion Bjj(t), t > 0, is Gaussian, has stationary increments, is self-similar with self-similarity 
exponent < H < 1, but is not Markovian, unless H = 1/2, in which case the fractional Brownian motion 
becomes Brownian motion. When 1/2 < H < 1, the increments of fractional Brownian motion have long-range 
dependence [49]. 

Because Brownian motion is Markovian with stationary increments, its finite-dimensional distributions can 
be obtained from the marginal density function 

= — Le"*"/ 4 ', xem (l) 

Vint 

at time t > 0. This density function is the fundamental solution of the "standard" diffusion equation: 

d t u(x, t) = d xx u(x, t), (2) 
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which in integral form reads: 



u(x, t) = Uo(x) + / d xx u(x,s)ds, uq(x) = u(x, 0). (3) 



(i 



Thus, Jb{x, t) is a solution of Eq. ([3]) with u (x) = S(x), where 5(x) is the Dirac delta distribution. We allow, 
throughout the paper, functions to be distributions. 

Remark 1.1. We follow the physics convention of not including the factor 1/2 in Eq. Q. Therefore, in this 
paper, "standard" Brownian motion B(t), t > 0, is such that, for each time t > 0, B(t) ~ iV(0,2i). The "tilde" 
notation X ~ fx (x) indicates that the random variable X has the probability density function fx (x) . 

Our goal is to extend Eq. ([3]) to non-Markovian settings. We will consider non-local, fractional and stretched 
modifications of the diffusion equation. These modified equations will be called Non-Markovian diffusion equa- 
tions, because, while they originate from a diffusion equation, the corresponding process, whose probability 
density function is a solution of these modified equations, will be typically non-Markovian. 

To motivate the modifications, consider first the non-random process l(t) = t, t > 0, which depicts a 
non-random linear time evolution and let fi(r,t) denote its density function at time t. Therefore one has 
fi(r,t) = S(t — t) where 6(x) is the Dirac distribution. It is natural to interpret fi{r,t) as the fundamental 
solution of the standard forward drift equation: 

d t u{T,t) = -d T u(T,t), T,t> 0, (4) 

which in integral form reads: 

u(r,t) = uq(t) — / d T u(T 7 s)ds 7 uq(t) = u(t, 0). (5) 
Jo 

The general solutions are of the form u(r,t) = uq(t — t) and thus, when uq(t) = S(t), the solution of Eq. 
is indeed u(t, t) = S(t — t). Observe that the variable r > plays the role of a space variable. 

We will consider the following generalization of the forward drift equation ((5|) 

u(r,t) = u (t) - [ K(t- s)d T u{T,s)ds, T,t>0, (6) 



where K(t), with t > 0, is a suitable kernel chosen such that the fundamental solution of Eq. J6]) is a probability 
density function at each t > 0. We refer to Eq. Q as the non-Markovian forward drift equation. 

The presence of the memory kernel K in Eq. ^ suggests a corresponding modification of the diffusion 
equation ([3]). Namely, we will consider the equation: 

u{x,t) = u (x) + [ K{t- s)d xx u(x,s)ds, t > 0. (7) 

Jo 

Its fundamental solution turns out to be: 

f( X ,t) = / G{x,T)h{T,t)dT, (8) 



where 

G(x.t) = -, 

/4nt 



G(x, t) = — = cxp(-x 2 /tt), (9) 



and h(r, t) is the fundamental solution of Eq. ([(j]). 
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The solution J8]) is a marginal (one-point) probability density function. We will consider different random 
processes whose marginal probability density function coincides with it. As illustration, consider the following 
example^]. 

Example 1.1. If we choose: 

K(t) = t -^, t>0, (10) 
V 71 " 

then we have, see Eq. J65]and Eq. ([70]) : 

Mr,t) = -^expf-^j , r>0, t>0, (11) 

as the fundamental solution of Eq. ([6]). Now consider the process 

D(t) = B(l(t)), i>0, (12) 

where B is a "standard" Brownian motion and l(t) > is a random time-change (not necessarily increasing), 
independent of B, whose marginal density function is given by h(r, t). One possible choice for the random time 
process is simply: 

l(t) = \b(t)\, t>0, 

where b(t), t > 0, is a "standard" Brownian motion [9,18]. Such a random time process l(t), t > 0, is self-similar 
of order H = l/2@ Let now B(t), t > 0, be another "standard" Brownian motion independent of b(t). Thus, 
the process (see also [3]) 

D(t) = B(\b(t)\), t>0, (13) 
has marginal density defined by Eq. ([Sj) with h(r,t) given by Eq. (fTTj) . 

But, D(t) is not the only process with density function f(x,t), given by Eq. ((8]). For example, the process 



Y(t) = VWPv 4 (t), t>0, (14) 

where -B1/4 is an independent fractional Brownian motion with self-similarity exponent H = 1/4, has the same 
one-dimensional probability density functions as the previous process D(t), i > 0, see Eq. (j40"l) with (5 = 1/2. 

Example 1.2. The fractional Brownian motion in Eq. (fl4)l has a self-similarity exponent H < 1/2. The 
increments of such a process are known to be negatively correlated [31,32,49]. To allow for the presence of 
fractional Brownian motion Bn(t) with < H < 1, we introduce a second (non-random) time-change t — ► g(t), 
where g(0) = and g(t) is smooth and increasing, that is we consider the non-Markovian diffusion equation 

u(x,t) = u (x) + / g'(s)K (g(t) - g(s)) d xx u(x, s)ds. (15) 
Jo 

whose fundamental solution is now: 

f(x,t)= I G(x,T)h(T,g(t))dT, (16) 
Jo 

where h is the fundamental solution of Eq. ([6]). If K(t) is as in Eq. (fT0)l and g(t) = t 2a , with < a < 2, then 
the processes: 

D(t)=B{\b{t^)\) 1 t>0, 



Y(t) = v 4Hrp a/2 (t), t>o, 

have a marginal density function defined by Eq. l[T6|) with h(r,t) as in Eq. (fTTj) . which is the fundamental 
solution of Eq. (fTB"]) . In this case V(t) is defined through an independent fractional Brownian motion B a / 2 with 
Hurst's parameter H = a/2 and thus < H < 1. This is a special case of Eq. l[78"l) . 

1 In these examples we refer to facts which are justified later in the paper through forward references. The reader may want to 
focus at this point only on the examples and ignore the references. 

2 Another possible choice for a random time process with marginal density given by Eq. is the local time in zero of a 

"standard" Brownian motion [4]. In this case the time-change process l(t) is increasing. 
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The preceding examples illustrate the themes pursued in the paper. We will focus, however, not only on 
power- like kernels such as those defined in Eq. lfT0|) . but also on exponential- like kernels such as: 

K(t) = e- at , a > 0. (17) 

We also consider what happens when the Brownian motion B(t), t > 0, is replaced by a more general linear 
(time-homogeneous) diffusion Q(t), t > 0, governed by the Fokker-Planck equatiorll, 

d t u{x,t)=P x u(x,t), (18) 

where V x is a linear operator independent of t acting on the variable x G K. In other words we consider the 
non-Markovian diffusion equation: 

u(x,t) = u (x) + f g'(s)K(g(t)-g(s))P x u(x,s)ds. (19) 
Jo 

We show that its fundamental solution is: 

/■OO 

f(x,t) = / g(x,r)h(r,g(t))dT, (20) 
Jo 

where G(x,t) is the fundamental solution of Eq. (fT9| . while h(r,t) is the fundamental solution of Eq. |(6]). We 
also provide explicit solutions when V x is the differential operator associated with Brownian motion with drift, 
when it is associated with Geometric Brownian motion and when the kernel K(t) is the power kernel and the 
exponential kernel. 



In order not to dwell on technicalities, we suppose implicitly, throughout the paper, that we have sufficient 
regularity conditions, to justify the algebraic manipulations that are performed. The paper is organized as 
follows: 

• Historical notes are presented in Section [H 

• In Section [3] we study the non-Markovian forward drift equation ([6]) and its corresponding random time 
process l(t). We derive suitability conditions on the kernel K(t). We end the section by noting that a 
self-similar time-change process, for instance with self-similarity parameter H = (3, requires the choice 
K{t) = CtP~ 1 /T{f3) with < < 1. 

• In Section 2] we study the non-Markovian diffusion equation lfT5|) and its solutions, and we discuss its 
various stochastic interpretations. 

• In Section [5] we illustrate the fact that the stochastic representation is not unique. 

• In Section [6] we study the more general non-Markovian Fokker-Planck equation and derive its solution Eq. 
|20j). 

• In Section [7| we go thorough several examples with V x u{x,t) = d xx u{x,t), that is, when the underlying 
diffusion proces is Brownian motion. We consider non-Markovian diffusion equations, associated with the 
(3-power kernel K(t) — t /3_1 /r(/3), < (3 < 1, and with the exponential- decay kernel K(t) = e~ at , a > 0. 
We also consider different choices of the deterministic scaling function g(t), for example a logarithmic time 
scale g(t) = log(£ + 1) is considered. 

• In Section [8] we focus on applications when the underlying diffusion process is not standard Brownian 
motion. We consider the case of Brownian motion with drift and Geometric Brownian motion and we 
study the corresponding equations with the /3-power kernel and the exponential-decay kernel. 

• Section [9] contains a summary and concluding remarks. 

3 Also known as the forward Kolmogorov equation. 
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2 Historical notes 



Non-Markovian equations like Eq. ([7]) , or more generally Eq. (fT9|) , are often encountered when studying physical 
phenomena related to relaxation and diffusion problems in complex systems (see Srokowsky [47] for examples) . 

Equations of the type Q have been studied for example by Kolsrud [22]. He obtained Eq. [8]), but without 
providing specific examples. A similar study was done by Wyss [51] who, however, focused only on power-like 
kernels K(t) = Ct^ 1 . 

Sokolov [45] (see also Srokowsky [47]), studied the non-Markovian equation 

d t P(x,t)= [ k(t- s)L x P{x,s)ds, (21) 
Jo 

where L x is a linear operator acting on the variable x. He provided a formal solution in the form of Eq. |20|) . 
Observe, however, that our equation lfl9|) differs from Eq. (|2"Tjl . not only by the presence of the scaling function 
g(t), but also by the choice of the memory kernel. Our kernel K(t) and Sokolov's kernel k(t) are related by the 
equation: 

K{t) = I k{s)ds =S> K(s) = k(s)/s, s > 0, (22) 
Jo 

where the tilde indicates the Laplace transform, see Eq. lj25j) . The suitability conditions for these memory 
kernels are thus not the same (these conditions are developed in Section [3]). For example, consider the simple 
exponential- decay kernel e~ at , a > 0. This choice of the kernel is "safe" in the context of Eq. (fT9|) . i.e. for the 
choice K(t) = e~ at , but is "dangerous" if one considers Eq. (pT]) with the kernel k{t) = e~ at . In the case of 
Eq. (fT9|) . the exponential- decay kernel corresponds to a system for which non-local memory effects are initially 
negligible. In fact, Kit) = e~ at — ► 1 as t — > and thus the system appears Markovian at small times. On the 
other hand, the choice k(t) = e~ at corresponds to the kernel K{t) = a (1 — e~ at ) which for small times behaves 
like t. In this case Sokolov [45] noticed that the corresponding equations are only reasonable in a restricted 
domain of the model parameters and for certain initial and boundary conditions. 

Our starting point is different from that of the previous authors. Instead of starting directly from the Fokker- 
Planck equations ifl8)l . we start from the forward drift equation (0 which is then generalized by introducing a 
memory kernel K(t), Eq. ([6]). One is then naturally led to the non-Markovian diffusion equations lfT5l and (fT9l) 
after the introduction of the scaling function g(t). In fact, in specific cases, it is sometimes simpler to solve first 
the non-Markovian forward drift equation and then use the solution to solve the non-Markovian diffusion 
equation (fT5|) or (fT9|) by using ifTG)) or l(20|) . The form of the solution lfT6|) or ([20]) has now a ready-made interpre- 
tation. For example, in Eq. (f2Q|) the function Q(x, t) is the fundamental solution of the Markovian equation (fl8|) 
and the function /i(r, t) is the fundamental solution of the non-Markovian equation ([6]) and it is these two solu- 
tions that contribute to Eq. (|20|) which is the fundamental solution of the non-Markovian diffusion equation (fl9|) . 

Furthermore, the form ffl6|) or (J20j) has a natural interpretation in terms of subordinated processes, see Eq. 
(p"2"Tl . According to Whitmore and Lee [23], the term "subordination" was introduced by Bochner [5,6]. It refers 
to processes of the form Y(t) = X(l(t)), t > 0, where X(t), t > 0, is a Markov process and l(t), t > 0, is a 
(non-negative) random time process independent of X. The marginal distribution of the subordinated process 
is clearly: 

/>oo 

f Y (x,t) = / f x (x,r)fi(T,t)dT, t>0, xel, (23) 
Jo 

where fx(x,t) and //(r, t) represent the marginal density functions of the processes X and I. Therefore, Eq. 
(fTBl) or Eq. ([20)1 can be interpreted in terms of subordinated processes, with Eq. J6]) characterizing the random 
time process l{t) and Eq. ((THJ) characterizing the Markov parent process X(t). 
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The stochastic interpretation through subordinated processes, first suggested by Kolsrud, is very natural 
because Y(t) = X(l(t)) has a direct physical interpretation. For example, in equipment usage, X(t) can be the 
state of a machine at time t and l(t) the effective usage up to time t. In an econometric study, X(t) may be 
a model for the price of a stock at time t. If l{t) measures the total economic activity up to time t, the price 
of the stock at time t should not be described by X(t) but by the subordinated process Y(t) = X(l(t)). The 
resulting subordinated process Y(t) is in general non-Markovian. In this way, the non-local memory effects are 
attributable to the random time process l(t) and to its dynamics which is in general non-local in time, see Eq. ([6]). 

Note, however, that the solution of Eq. (fl9]l represents only the marginal (one-point) density function of the 
process and therefore cannot characterize the full stochastic structure of the process. As we note in the paper, 
there are also processes that are not subordinated processes that serve as stochastic models for non-Markovian 
diffusion equations like Eq. (fT9|) or Eq. (f2T]l . 

For example, consider in Eq. ([7]) the (3-power kernel K(t) = t 13 ^ 1 /T((3), with < (3 < 1. From a stochastic 
point of view, the fundamental solution of this equation, also called the time-fractional diffusion equation 
of order (3, can be interpreted as the marginal density function of a self-similar stochastic processes with 
parameter H = (3/2. This process, for example, can be taken to be a subordinated process Y(t) = B(l(t)), 
with a suitable choice of the random time I. In Kolsrud [22], the random time / is taken to be related to 
the local time of a d = 2(1 — /3)-dimensional fractional Bessel process, while in Meerschaert et al. [34] (see 
also [1,15-17,21,40,48]), in the context of a Continuous Time Random Walk (CTRW), it is chosen to be the 
inverse of the totally skewed strictly /3-stable process. The interested reader is referred to the wide literature 
concerning the relationship between CTRW and non-Markovian diffusion equations and its applications. See 
for instance, [2,13,14,19,28,35,36,41,42,50,52] and references therein. 

Schneider [43], moreover, in a very general mathematical construction, introduced the so-called Grey Brow- 
nian motion. This process is a self-similar process with stationary increments which, as turns out, can be 
represented by Y(t) = ApBn(t), t > 0, where Bh is a fractional Brownian motion with H = f3/2 and is a 
suitable chosen random variable independent of Bh (see Mura et al. for details [38,39]). This process has a 
marginal density function that evolves in time according to the time-fractional diffusion equation of order (3. In 
this case the non-Markovian property is due to the presence of the fractional Brownian motion. As we show in 
the paper, long-range dependence can be made to appear through the time-scaling function g(t), see Eq. (fl"5| 
and Example 11.21 Figures 4, 5 and 6 display trajectories of the processes D(t) and Y(t) and corresponding 
density functions. 

3 The non-Markovian forward drift equation 

We start with the following generalization of Eq. |(5|), namely: 



where K(t), with t > 0, is a suitable chosen kernel. We then choose a random time process l(t) such that, 
for each t > 0, its marginal density fi(r,t) is the fundamental solution of Eq. |24|) . Observe that Eq. |24|) is 
"non-local" because u(r, t) involves u(t, s) at all < s < t. Equation ([24|) will be called non-Markovian forward 
drift equation, see Section [TJ Eq. ([6]). 

It is convenient to work with Laplace transforms. We indicate by L{tp(x, t); t, s} the Laplace transform of 
the function <p with respect to t evaluated in s > 0, namely: 



If the function ip depends only on the variable t we write simply ip(s), because in this case there is no ambiguity 
concerning the integration variable. In particular we let K (s) denote the Laplace transform of the kernel K. 




(24) 




(25) 



£{/,(M);i,s} = -J— expf --J— V r,s>0, (26) 



Proposition 3.1. Let fi(r,t) denote the fundamental solution of Eq. \21$ . Then 

1 / r 

— — exp - — — 

sK(s) \ K(s) 

and zero for r < 0. 

Proof: we take the Laplace transform with respect to the variable t in Eq. (|24|) : 



9t £ (T)S) = ^_3i4 (27) 

thus Eq. ([26]l is a solution, in the distributional sense, when u (t) = S(t). Indeed the general solution ofEq. 
l27j) with u q (t) = 5(t) is: 

ip(r,s)= ~ exp ( -^J— I +Cexp ( -^J— V rel, 
where C is a real constant and where 

= { j; :%° ' (28) 

is the Heaviside's step function. Since we require ip(r, t) = for r < 0, we get C = i.e. Eq. (J26J) . □ 



3.1 Suitability conditions on the kernel K 

We must choose the kernel K such that the fundamental solution of Eq. ([2%)l is a probability density in r > 0. 
We observe that if fi(r,t) satisfies Eq. lj24|l and Eq. (|26j) . then it is automatically normalized for each t > 0. 
In fact, for a function (^(x, t) for which it is always possible to change the order of integration, one has: 

<p(x, t)dx = 1 <==> / tp(x, s)dx = s~ . (29) 



Since Eq. ([261 satisfies the right-hand side of Eq. (|29|). we get J R+ //(t, t)dr = 1. One still needs, however, to 
choose the kernel K such that /;(t, i) > for all r, t > 0. 

In order to get a suitable condition on the kernel K, we make use of the notion of completely monotone 
function. Recall that a function (f(t) is completely monotone if it is non-negative and possesses derivatives of 
any order and: 

(~l) k — p(t) >0, *>0, fceZ + = {0,1,2,...}. (30) 

We observe that as t — > 0, the limit of d k (p(t)/dt k may be finite or infinite. Typical non-trivial examples are 
<p(t) = cxp(— at), with a > 0, ^(t) = 1/i and <fr(t) = 1/(1 + t). It is easy to show that if ip and -0 are completely 
monotone then their product (pip is as well. Moreover, if ip is completely monotone and ip is positive with first 
derivative completely monotone then the function y(^) is completely monotone. 

We have the following characterization of completely monotone functions [10]: 

Lemma 3.1. A function (p(s), defined on the positive real line, is completely monotone if and only if is of the 
form: 

poo 

<p(s) = / e- ts F{dt), s > 0, 
Jo 

where F is a finite or infinite non-negative measure on the positive real semi-axis. 

Hence, to ensure that fi{r,t) > for all r, t > 0, it is enough to require that the function defined in Eq. 



([26]) must be completely monotone, as a function of s, for any r > 0, and thus that the kernel K satisfies the 
following: 
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Suitability conditions 

1. sK(s) is positive with first derivative completely monotone, 

2. 1/K(s) is positive with first derivative completely monotone. 

Indeed, we can view Eq. [26)1 as the product of the two completely monotone functions 1/u and exp(— tu), the 
first evaluated at u = sK(s) and the second evaluated at u = 1/K(s). 

3.2 Examples 

Example 3.1 (fi-power kernel). If we choose: 

m " WY 

we get K(s) = s" 13 . In this case sK{s) = s 1 ^ 13 is positive and has first derivative (1 — f3)s~ 13 completely 
monotone if and only if < /3 < 1. Moreover, 1/K(s) = s 13 is positive with first derivative /3s^ _1 completely 
monotone if and only if < f3 < 1. Therefore, a good choice for the kernel K is: 

K ( t ) = Wfi> 0</3<l. (31) 
Example 3.2 (Exponential- decay kernel). Choosing: 

K{t) = exp(-ot), a>0, (32) 

we get sK(s) = s/(s + a) which is positive with first derivative a(s + a)~ 2 completely monotone for any a > 0. 
Moreover, 1/K(s) = (s + a) is positive if a > with first derivative completely monotone. 

Example 3.3 (/3-power with exponential- decay kernel). Choosing: 

K{t) = ^— exp(-at), < < 1, a > 0, (33) 

we have K(s) = (s + a) _/3 . Therefore, sK(s) = s(s + a)~P which is positive if a > with first derivative 
(s + a)~^(l — (3s(s + a) -1 ) completely monotone if < /3 < 1. Moreover, 1/K(s) = (s + a) 13 is positive if a > 
with first derivative f3(s + a) 13 ^ 1 completely monotone if < (5 < 1. 

The following theorem states that a self-similar random time process l(t), t > 0, is associated with the kernel 
K(t) in Example O 

Theorem 3.1. // the time-change process l(t), t > 0, is self-similar (for instance of order H = (3), with 
marginal probability density fi(r,t) satisfying Eq. h26\) . then we must have: 

K[t) = G wr 0</3 - 1 ' (34) 

for some positive constant C . 

Proof: The self-similarity condition entails that for any t, t > and for any a > 0: 

a- f3 f l {a- r,t) = fi{r,at). 
If we take the Laplace transform and set /(r, s) = C{fi(r, t);t, s}, we have: 

a-^(a-"r,«) = ^(r^). 

a V a/ 



8 



Using Eq. l26j) we get that for any r, s > and a > 0: 

a-' 3 / ffl ~M 1 / t_A 

Since this relation is valid for any choice of r > and s > 0, putting r = and s = a, we get: 

p-0 1 
#(a) ~~ 

Thus, for any a > 0: 

if (a) = £(l)a _/J , 

which is the Laplace transform ofEq. l34|) . If we add moreover the condition of complete monotonicity we find: 
< (3 < 1 as indicated in Example 13,11 □ 



4 Non-Markovian diffusion equation 



We focus here on the non-Markovian diffusion equation (|T5j) introduced in the first section. There are two 
ingredients: 

1. The fundamental solution ofEq. |24| . denoted here by h{r,t) and defined by Eq. ((26j) . 

2. The fundamental solution G(x,t), defined by Eq. J9]), of the standard diffusion equation which is the 
one-dimensional density of the "standard" Brownian motion. 

The following theorem combines these two ingredients and provides the fundamental solution of a corresponding 
non-Markovian diffusion equation. 

Theorem 4.1. Let h(r,t) denote the fundamental solution of Eq. \2$ , so that by Proposition ^. 1\ one has: 

L{h(T,t);t,s} = — cxpf-^— I, t, s>0, (35) 
sK(s) \ K(s) J 

for a suitable choice of K. Let g be a strictly increasing function with g(0) = and let G(x,t) be defined by Eq. 
J3). Then, 



f(x,t)= / G(x,r)h(r,g(t))dT, (36) 
Jo 

is the fundamental solution of the non-Markovian diffusion equation: 

u{x, t) = u (t) + I g'(s)K (g(t) - g(s)) d xx u{x : s)ds. (37) 
Ja 

Proof: see Section El □ 

We have immediately the following: 

Corollary 4.1. Lf H(x,t) is a solution of the standard diffusion equation with initial condition H(x,0) = uq(x), 
then the function: 

/>oo 

u(x,t)= H{x,T)h(T,g(t))dr (38) 







is a solution of Eq. $37)) . 
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Proof: If, for any t > 0, the function f(x, t) defined in Eq. ([36)) is the fundamental solution of Eq. <f37|l then a 
general solution is given by: 



u(x,t) = / f(x - y,t)u (y)dy = / / G(x - y,T)u {y)h{r, g(t)) dr dy 
Js. Jr jo 

G(x-y,r)u (y)dy]h(T,g(t))dT= [ H(x,T)h(r,g(t))dT. □ 



/ Jo 
We observe that: 

1. The equation f35|) states that h(r,t) is the fundamental solution of Eq. ([24"1) . 

2. While G(x,t) is the fundamental solution of the standard diffusion equation obtained when u (x) = S(x), 
the general solution, denoted H (x, t) in the above theorem, results from a general initial condition uq{x). 

Many physical phenomena, especially related to relaxation processes in complex systems, are described by non- 
Markovian "master equations" like Eq. <f37|l . K (t) is a memory kernel and g(t) is just a "time-scaling" function. 
Such equations are often argued by phenomenological considerations and can be more or less rigorously derived 
starting from a microscopic description [7,20,47,53]. 

5 The stochastic representation is not unique 

The solution of the non-Markovian diffusion equation can be viewed as the marginal density function of the 
subordinated process, see Eq. 11121) 

D(t) = B(l(g(t))), t>0, 

since its marginal density is: 

/"OG 

f D (x,t) = / G{x,T)fi{T,g{t))dr. 







Here, for each t > 0, D(t) ~ f D (x,t), B(t) ~ G(x,t) and l(t) ~ fi(r,t). In the notation of Theorem 14. 1[ we 
have fD{x,t) = f(x,t) and fi(r,t) = h(r 7 t). The Laplace transform of /;(t, t) with respect to t is given by Eq. 

131. 



This stochastic representation is not unique (see Example 11,14 Example 11.21 and examples below). Indeed, 
the non-Markovian diffusion equation characterizes only the marginal, that is one-point, probability density 
function. However, processes with a different dependence structure can have the same marginal density f(x,t). 
Additional requirements could be imposed so as to specify the stochastic model more precisely. 

Example 5.1. If we require the random time process lp{t), t > 0, to be self-similar of order (3, then in view of 
Theorem 13. 1\ the kernel must be chosen as in Eq. |34^ and we must have < (3 < 1. We will study this case 
more in details in Section [71 Here we just observe that if we consider a "standard" fractional Brownian motion 
of order /3/2, then f{x,t) is also the marginal distribution of 



Y(t) = yJl^X)B m {t), (39) 
where Bp^it) is assumed to be independent of Z/3(l). 
In fact, because lp(t), t > 0, is self-similar of order H = 0, one has: 

D(t) - B(l p (t)) = y/wfiBO) = y^(l)^/ 2 S(l) = ^T)0 2 B m (l) = yJl^B^t) - Y(t), (40) 

where = denotes here the equality of the marginal distributions. 

Both D(t), t > 0, and Y(t), t>0, are self-similar processes with Hurst's exponent H = (3/2. However, while 
Y(t), t > 0, has always stationary increments, this is not in general true in the case of the process D(t), t > 0. 
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6 Non-Markovian Fokker-Planck equation 

We considered up until now processes of the type B(l(g(t))), where B is a "standard" Brownian motion. What 
happens if we replace B by a more general diffusion? Namely, what happens if instead of starting with the 
standard diffusion equation ((2]) we start with a more general Markovian Fokker-Planck equation: 

d t u(x,t) = V x u(x,t), t>0, (41) 

where V x is a linear operator, independent of t, acting on the variable xl We have the following generalization 
of Theorem 14.11 

Theorem 6.1. Suppose that h(r,t) is a probability density function satisfyingEq. [2B\) 

Z{h(r,t) i t,8} = -^exp(--J—), r, S >0, (42) 
sK(s) \ K(s) J 

for a suitable choice of K . Let g be a strictly increasing function with g(0) = and G(x, t) be the fundamental 
solution of Eq. H41\ ). Then the fundamental solution of the integral equation: 

u(x,t)=u (t)+ [ g'{s)K(g{t)-g( S ))V x u{x,s)ds (43) 
Jo 

is 

f{x,t)= / G(x,r)h(T,g(t))dT. (44) 
Jo 

We provide two versions of the proof. The first starts with the solution f(x,t) in Eq. (|l4f and verifies that it 
satisfies Eq. (03]). The second starts from the partial integro-differential equation l[4*3|) and derives the solution 
f(x,t) under certain assumptions stated below Eq. |49|). 

Proof 1: For first we observe that 

L{f(x,t);g(t),s} = -J—LiGix^^t^Kis)- 1 }. (45) 
sK(s) 

With the change of variables g(s) = z, we write: 

u(x, g~ (w)) = Uq(x) + / K (w — z) V x u(x, g~ 1 (z))dz, w = g(t). (46) 
Jo 

We want to show that Eq. I|44p with the choice (|42| solves Eq. (|43f . If we take the Laplace transform of Eq. 
(03J) using Eq. we get: 

L{u(x,t);g(t),s} = + K(s)V x L{u(x,t); g(t), s} 

s 

that is: 

s£{u(x, t);g(t), s} ~ u {x) = sK(s)P x L{u(x, t);g(t), s}. (47) 
Now, if we substitute on Eq. (|47f a solution of the form ((44]) . 



u{x,t)= H(x,T)h(r,g(t))dT, (48) 
Jo 

we have: 

Kis^Linix, t); t, K(s)- 1 } = u {x) + V x L{H{x, t);t, Kis)^ 1 } 
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i.e. we have, with obvious notations: 

tH(x,t) — u (x) + V x H(x,t), 

in which one readily recognizes the Laplace transform of the Markovian Fokker-Planck equation with the same 
initial condition uq(x). Therefore: 

d t H(x,t) = V x H(x,t), H(x,0) = u (x). 



This argument shows not only that Eq. l(44|) is the fundamental solution of Eq. l(43j) . but also that a general 
solution is given by Eq. l(i8|) (see Corollary 14. ip . This result is summarized in Corollary 16.11 (see below) . 

Proof 2: We now start from Eq. lj43j) and we use integral transforms in order to get the fundamental solu- 
tion. Let T denote the Fourier transform operator and let: 

Akx , 



{Tip){k,t) = (p{k,t) = / e lkx ip{x,t)dx. 

Since Uo{k) = 1, and since {TT> x u)(k,t) = {TV x T- x Tu)(k,t) = V k u(k,t), where V k = {TV X T- X ) k denotes the 
Fourier transform of the operator V x , we have: 

u{k,g^ 1 (w)) = 1 + / K(w - z)T k u(k,g- 1 (z))dz. 



Taking the Laplace transform we have: 

H{u(k, g~ 1 (w));w, s} = s^ 1 + V k K(s)L{u(k, g^ 1 (w));w, s}, 

which is the same as: 

H{u(k, t);g(t), s} = s- 1 + V k K(s)L{u(k, t);g(t),s}. 

Therefore: 

(K(s)- 1 - V k ) H{u(k, t); g(t), s} = s" 1 ^)" 1 . 

Denoting l(fc) = 1, we have: 



C{u(k,t);g(t),s}=^— (kis)- 1 -^) \(fc), (49) 
sK(s) v ' 

where we suppose that the operator (^K(s)^ 1 — Vkj is well defined and acts on the constant function l(fc) = 1. 

Observe that the Fokker-Planck equation ((41]) is obtained from Eq. lj43"|) by setting K(t) = 1, for each t > 0, 
that is K(s) = s , and g(t) = t, for each t > 0. In this case Eq. (f49| becomes: 

L{g(k,t);t,s} = (s -Vk)- 1 !^). (50) 
where Q{x,t) is the fundamental solution. Taking the inverse Fourier transform, we get: 

L{g(x,t);t,s} =J- 1 {(s-V k )- 1 l(k) ; k,x}, (51) 

where: 

^{ipfcs) ; k,x} = -!- / e- ikx ip(k,s)dk. (52) 
Replacing s by K(s)~ 1 in Eq. l(51"j) . one has: 

L{g(x,t);t,K(s)- 1 } = ?- 1 {(K(s)- 1 -V k )- 1 l(k); k,x} . (53) 
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Going back to Eq. (|49| and inverting the Fourier transform we obtain in view of Eq. ([53]) : 

L{u(x,t);g(t), S } = -^—^ ( (k(s)- 1 - V k ) ~* l(fc) ; k,x\ = -J— L{G( X> *); i, ^(s)" 1 }. 
sK(s) I v 7 J sA (s) 

that is Eq. (45]). □ 



5 2 

Remark 6.1. If the Markovian process is a Brownian motion one has V x = tt~o- The Fourier transform of V x 

^ ox 2 

is Vk — -k 2 and Eq. ((49J) becomes: 



H{u(k,t);g(t),s} = -J— (Kisy'+k 2 ) 1 l(fc), 
sK(s) v ' 



where 

-l 



(K(s)-' + k 2 ) i(k) = TZ T 

V > (K(s)- 1 + fc 2 J 



which is well defined because K(s) 1 is positive. 

Corollary 6.1. If H(x,t) is a general solution of the Markovian Fokker-Planck equation H41\ ) with initial 
condition 7i(x,0) = uq(x), then the function: 

r ^ 

u(x,t)= H(x,T)h(T,g(t))dT (54) 



is a general solution of Eq. j43\) . 

>From a stochastic point of view, f(x, t) could be seen as the marginal distribution at time t of the subor- 
dinated process: 

V(t) - Q(l(g(t))) (55) 

where Q is the diffusion governed by the Fokker-Planck equation f4T]) and l(t) is the random time process, 
independent of Q(t), with marginal distributions defined by h(r,t). 



7 Examples involving standard Brownian motion 

In the following examples, we consider stochastic models where the operator V x in Eq. (|4*T]) is d xx , namely 
the operator corresponding to standard Brownian motion. We will study more general operators in the next 
section. We shall choose various kernels K(t) and various stretching functions g(t). We let h(r,t) denote the 
fundamental solution of the non-Markovian forward drift equation (]24| . Since the corresponding stochastic 
models are not unique, we will mainly focus on the subordinated process B(l(t)), t>0. However, we also give 
examples of other appropriate stochastic models. 



7.1 Time-fractional diffusion equation 

Let g(t) = t. Consider the (i-power kernel: 

K{t) = Wy 0<p - h (56) 

and let h(r, t) denote the fundamental solution of the non-Markovian forward drift equation (|24]l with kernel 
ifHfil. 



Remark 7.1. In view of Theorem 13.11 such a kernel arises if one requires h(r, t) to be the marginal density 
function of a self-similar random time process l(t) of order /?. 
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InsertingEq. l56j) in Eq. lj37j) we obtain the following equation: 



u(x, t) = Uo(t) + 



(t - sf 1 d xx u(x, s) ds . 



(57) 



which is sometimes called the time-fractional diffusion equation [27,44]. In view of Theorem l4.14 the fundamental 
solution is: 



f(x,t) = / G(x,r)h(r, t) dr, 



where h(r, t) satisfies: 
Such a function h(r, t) can be expressed as: 



o 

Z{h{T,t);t,s} = s l3 ~ 1 e- TS ' 



t,s > 0. 



h(T,t) = t- Mp(Tt-P), 
where Mp(r), is defined for < (3 < 1 by the power series [24,25]: 

(-r) k 



(58) 
(59) 



k=0 1 



-f3k + (1 - /?)] 



^^^rp(fc + l))]sin[^/3(fc + l)], 



(60) 



7T A — ' k\ 

fe=0 



The above series defines a transcendental function (entire of order 1/(1 — /?)) [12]. 




Figure 1: Plot of the density function h(r,t) = t ®M{rt ^) at time t = l, for different values of the parameter 
0= [1/4,1/2,3/4]. 
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Remark 7.2. The function h(r, t) in Eq. (|59| represents the fundamental solution of the time-fractional 
forward drift equation (see also [15]): 

1 f f 

u(t, t) = uo(r) - p^y ^ (* - sf- x d T u{r, s)ds. (61) 
This equation reduces to the standard drift equation when (3 — ► 1. 
7.1.1 Properties of the M-function 

It is useful to recall some important properties of the M-function [12,29]. These are best expressed in terms of 
the function 

Mp(T,t)=t-PM f3 (Tt-P), (62) 

defined for any r, t > and < (3 < 1. 

1. The Laplace transform of M.p(r, t) with respect to t is: 

H{Mp(T,t);t,s} = s^V™", r, s > 0. (63) 

2. The above equation suggests that in the singular limit j3 — > 1 one has: 

Wi(T,t) =(5(r-t), r,t>0. (64) 

3. K/3= 1/2: 

M 1/2 (T,t) = -J=cx P (-t 2 /4*), r,t>0. (65) 
V 7rf 

4. The M-function is a particular case of a Fox i/-function [30,44]. We indicate with 

/>oo 

M{ip(x);x,u} = / ^a^dx, (66) 
Jo 

the Mellin transform of a function ip(x), x > 0, with respect to x evaluated in u > 0. The Fox if -function 

TTTn.n/ \ rrm,n ( { a i 7 i=l , . . . ,p \ 

M M V (&i,&)i=l,...,, )' 

is characterized by its Mellin transform as follows: 



li ™ [z ^ z ^ C(u)D(u)- 



(67) 



with 



A(u) = JJr(6,- + &•«), B(u) = JJr(l 0j - «,•«), 

C(it)= r(l-&j-/3ju), D(u)= ' r( aj +aju). 

i— m+1 j— n+1 

Here: 1 < m < q, < n < p, aj,f3j > and a,j,bj G C (see [11,33,46] for more details). 



Starting from Eq. f63|) and skipping to the Mellin transform, it is easy to show that we have the following 
relation: 

M p (t, t) = t-f>H$ (rt-P\ ' 3) ) , r,t > 0, < (3 < 1. (68) 
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5. Using the representation Eq. ([68]) and Eq. ([67)) we have for any 77, /3 € (0, 1), see also [29]: 

poo 

M„(x,t)= / M v (x,T)Mp(T,t)dT, v = r)/3 x>0. (69) 



The expression (|59| for the function h(r,t) follows from Eq. {63]), that is: 

h(r,t) = Mp(r,t), r,t>0. (70) 

Moreover, when (5 — > l,Eq. [64]) gives h(r,t) = 5(t — t) as expected (see Remark l7,2|) . Comparing Eq. ([9]) and 
Eq. ([651 one observes that: 

G(x,t) = ±M 1/2 (\x\,t). (71) 

Using Theorem 14.11 and Eq. ([69]) together with Eq. ([?0]l and Eq. ((71)1 we recover the fundamental solution of 
the time-fractional diffusion equation [27]: 

f(x,t)= G(x,T)h(r,t)dT = - M 1/2 (\x\ 7 T)M p (T 7 t)dT 
Jo 1 Jo 

= \M m {\x\,t) = \r^Hl N2 {\x\t^l^). (72) 

Several plots of the M-function are presented: in Figure[T]the function h(r,t) = Mp(r,t) is drawn at a fixed 
time t = 1 and for different values of the parameter /?; in Figure[2]is presented the plot of f(x,t) = i^Vl /3/2 (I ^ I j 
at a fixed time t = 1 and for different values of (3\ in Figure [3] is shown the time evolution of f{x,t) for fixed 
P = 1/2. 




Figure 2: Plot of the density function f{x,t) given by Eq. fT2]l at time t = 1, for different values of the 
parameter [3 = [1/4, 1/2, 3/4, 1]. For (3 — 1 one recovers the standard Gaussian density ([7T]) . 
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Figure 3: Plot of the density function f(x, t) for fixed (3 = 1/2, at different times t = [0.1, 1, 10, 10 2 ]. 
7.1.2 Stochastic interpretations of the solution 

From a stochastic point of view, the function h(r, t) in Eq. I|59p can be regarded as the marginal distribution of 

hit), *>0, 

where lp(t), t > 0, is an H-ss random time with H = (3. We have that for each integer m > 0: 
In fact, from Eq. (|58|) . for each integer m > 0, we have : 



oo 



m s^- 1 e- TS ' 3 'dr = mis-" - 1 , 



T 







which, inverting the Laplace transform, gives Eq. (|73|) . 

For instance, with the suitable conventions [8], lp(t), t > 0, can be viewed as the local time in zero at time t 
of a d = 2(1 — /3)-dimensional Bessel process [37]. The function f(x,t) in Eq. lj7"2]) is then the marginal density 
function of 

D(t)=B(l (t)), 

which is self-similar with H = (3/2. In this case, because lp(t) is self-similar of order f3, we immediately have 
an example of a different process with the same marginal distribution of D(t) (see Example 15 .ip . In fact, if 
we consider a "standard" fractional Brownian motion Bp/ 2 of order (3/2, then f(x,t) can also be seen as the 
marginal distribution of 

Y(t) = ^)B 0/2 {t), (74) 
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Figure 4: Trajectory of the process B(lp(t j) (top panel), with < t < 1 and /3 = 1/2. The random time process 
is chosen to be h/2{t) = \b(t)\ where b(t) is a "standard" Brownian motion (see Example 1 1.1|) . The corresponding 
trajectory of the random time process is presented in the middle panel. The estimated variance, computed on 
a sample of dimension TV = 5000, is presented in logarithmic scale in the bottom panel and fits perfectly the 
theoretical curve 2i 1 / 2 /T(3/2). 

where Bpi 2 (t) is assumed to be independent of (see Example 15. ip . The process Y(t), t > 0, is called grey 
Brownian motion [43]. 



From Eq. (|73| one can derive immediately all the moments for the processes D(t) and Y(t). For any integer 
m > 

' E{D{t) 2m+1 ) = E(Y(t) 2m+1 ) = 0; 

2m' o ( 75 ) 
E(D(t)*») = E(Y(t) 2m ) = — ; 
^ r^m + lj 

Because < /3 < 1, the variance grows slower than linearly with respect to time. In this case one speaks about 
slow anomalous diffusion. Moreover, the increments of the fractional Brownian motion Bp/ 2 (t) do not have 
long-range dependence. In contrast, the next example allows for the presence of long-range dependence through 
the introduction of a scaling function g(t) = t a '^ (see also Example 1 1.2ft . 

7.2 "Stretched" time-fractional diffusion equation 

If in the setup of Section [7jTJ where the kernel K(t) is given by Eq. (|56|l . we introduce a scaling time 

g{t) = t a 'P 
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Figure 5: Marginal density function f(x,t) = ^A4i/i(\x\,t) of the process B{li/ 2 (t)) at time t = 1 and x G 
[—5,5]. The histogram is evaluated over N = 10 4 simulated trajectories of the process B(\b(t)\) (Figure[4|). 

with a > 0, then the integral equation ([571 is replaced by Eq. (|37|) . namely 

u(x,t)=u (t) + —- s0 1 (t0 -s?) d xx u{x,s)ds. (76) 
Therefore, using Eq. l59|) : 

ft(r, fl (t)) - g(t)-^Mp(rg(t)-^) = t- a Mp(rt- a ), 
and, using Eq. lf72"T) . the fundamental solution f(x,t) of Eq. fT6|) reads: 

/(x, t) = f(x, g(t)) = ±t- a / 2 M m (\x\t- a / 2 ), t > 0. (77) 

The function f(x, t), t > 0, is the marginal distribution of the process 

D{t) = B (h(t a ^))) , t > 0. 

The time-change process l/3(t a ^) is self-similar of order H = a and the process -D(i) is then self-similar with 
H = a/2. In the case < a < 2, the function f(x,t) is also the marginal density of 

y(t) = \fW) B c*/2( t ), *>0, 0<a<2, (78) 

where B a / 2 {t) is a "standard" fBm of order -ff = a/2 independent of The process 3^(i), t > 0, is called 

generalized grey Brownian motion [38]. 
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f(x,t) 




4 x 5 



Figure 6: Trajectory of the process y/lp{l)B a / 2 {t) (top panel), with 0<i<l,/3=l/2 and a = 3/2. The 
random variable £1/2(1) is Gaussian, see Eq. lj65"T) . The estimated variance, computed on a sample of dimension 
N = 5000, is presented in logarithmic scale in the middle panel together with the theoretical curve 2i 3 / 2 /r(3/2). 
In the bottom panel the histogram, evaluated over a sample of N = 10 4 trajectories, fits the exact marginal 
density Eq. ([77| at time t = 1. 



In this case, for any integer m > 0: 

E(D(t) 2m+1 ) = E(y(t) 2m+1 ) = 0; 

2to! 



E(D(t) 



2m\ 



E(y(t) 2m ) 



(79) 



r(/3m+l) 



We have slow diffusion when < a < 1 (the variance grows slower than linearly in time) and fast diffusion 
when 1 < a < 2 (the variance grows faster than linearly in time) . In this case the increments of the process 
y(t) exhibit long-range dependence. 



7.3 Exponential-decay kernel 

Let g(t) = t. With the exponential-decay kernel: 

K(t) = exp(-at), a > 0, t > 0, 

we obtain the following equation: 

u(x,t) = Uo(x) + / e _a(t_s ''9 2 ; a ;u(a;, s)ds. 



(80) 



(81) 



In this case K(s) = (s + a) 1 and the marginal distribution of the random time process l(t), t > 0, is defined 
by Eq. If26j) : 



( C{/ I (T,i);t,«} = i±^e-^ + °) I 
s 



T > 0. 
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Figure 7: Plots of the marginal density of the random time /j(t, t)Eq. [82]) as a function of r at times t = 
[0.5, 1, 1.5], and with a = 1. The vertical line corresponds to a point mass (delta function). 



Therefore, 

fl(r,t) = er Ta {5{t - t) + ad{t - t)) = e- ta S(r - t) + ae- Ta 6(t - r), (82) 

where 9{x) is the step function (|28|) . A graphical representation of the time evolution of fi(r,t) is presented in 
Figure [7J 

Remark 7.3. The function //(t, t) defined in Eq. [82| is the fundamental solution, in the sense of distributions, 
of the "exponential" forward drift equation: 

u(t, t) = u (t) - f e- a{t - s) d T u{T, s)ds. 
Jo 

This follows from Proposition 13. II To check it directly we note that fi(r,0) = S(t) and for any t > 0: 

t pt 
e- a{t ~ s) d T fi{T, s)ds = - e- a(t ~ s) d T (e- as S( T - s) + ae- aT 6(s - r)) ds 
Jo 

= - f e - a( *- s) (-e- as 5'{T -s) + ae- aT 6(s - r) - a 2 e - aT 9(s - r)) ds 
Jo 

= e~ at 5(T -t) + ae- at 9(t - r) + ae~ a{t+T ^e{t - r)(e at - e aT ) 
= e- at S(r -t) + ae- aT 9(t - r) = /,( T) t), 
where we have used the fact that: 

t 

5'{T-s)ds = 5(t-r). 

We observe that when o^Owe recover the forward drift equation ((U and indeed /z(r, t) = 6(t — t). 
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As noted in Example 13. 2\ Eq. ((82j) actually defines a probability density for any t > 0. The following proposition 
provides its moments. 

Proposition 7.1. For each integer m > one has: 

E(l(t) m ) = — (1 - e- at ) + e- at ft"* - ^ _Lt*a* _m J ■ ( 83 ) 
Proof: for any t > 0, we must evaluate: 

r m / ; (r, t)dT = e- at t m + a 

where we have used Eq. i|82|) . In order to evaluate the exponential integral in the above equation we write: 

a fr-e-^dr = (-l) m ad™ [(1 - e^){a^)} = (-l) m af^ (™) d k a (l - e~ at )d^ (a" 1 ) 
Jo k=0 ^ ' 

k=0 k=l 

thus one has Eq. EH}. □ 

The function //(r, t) can be written: 

/,(T,t) = e- at ,5(r-i) + (l-e- at )^(T,t), r,i>0, a > 0, (84) 

where: 

<p(r,t)=a 1 _\_ at , r,t>0, a>0. (85) 

Because fi(r,t) is a probability density, then so is (p(r,t). The corresponding random time process t > 0, 
can then be chosen to be: 

l(t) = b t t+(l-b t )j(t), *>0, (86) 

where &t, i > 0, is a stochastic process such that, for any fixed t > 0, bt is a Bernoulli random variable with 
Pr(b t = 1) = e _a * and Pr(bt = 0) = 1 — e~ at , and t > 0, is a stochastic process, independent of bt, with 
marginal distribution given by ip(r,t). 

Remark 7.4. The random time l(t) defined by Eq. (f86|) cannot be increasing everywhere. This is due to the 
fact that b t and j(t) are independent and Pr(j(t) < t) = 1 for any i > 0. Indeed, suppose that is increasing. 
This implies that for any t > and e > 0: 

1 = Pr(l(t + e)> l(t)\ b t = l) = Pr(l(t + e) > t) 

= Pr(l(t + e)>t\ b t+e = l)Pr(b t +e = 1) + Pr(i(i + e) > i| 6 t+e = 0)Pr(&t +£ = 0) 

= e -a(t+«) + A _ e -a(t+e)\ p^y. + e ) > f) 

= 1 - (l - e -°( t+e >) Pr(j(* + e) < t) 
therefore taking e^Owe get 1 = e~ at with a, t > 0, which is a contradiction as soon asa^O and t > 0. 

On the other hand, a trivial example of an increasing process with marginal distribution given by Eq. ((82)1 

is: 

7(t) = min(X,t), t > 0, (87) 
where X is an exponentially distributed random variable: X ~ ae~ aT , r > 0. 

We now turn to Eq. [ST]) . We have the following result: 
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Proposition 7.2. The fundamental solution of Eq. i81)) is: 

f(x, t) = e~ at G{x, t) + (1 - e~ at )(j)(x, t), (88) 

with: 

fa, *) = 4(1 ^_ at) {e^Erf f + - e -^Erf ^ - ^ - 2 sinh(|x|^)| , (89) 
where Erf(x) = -3= dy and where Erf(— x) = — Erf(ai). 



f(x,t) 0.8 




X 



Figure 8: Plot of the fundamental solution f(x,t), Eq. ([88]). at time t = 1, for different values of the parameter 
a = [0, 0.1, 1, 2]. When a = 0we have the standard Gaussian density. 

Proof: by Theorem 14.11 andEq. l85|) . the fundamental solution of Eq. ((8Tj) is: 

/>00 

f(x,t)= G(x,T)f l (T 7 t)dT = e~ at G(x,t) + (l-e- at )cj ) (x,t) 7 (90) 



where 

We have that: 

^fe *1 = 

1 - e~ at 



*) = / r)e- aT 9(t - r)dr. 



One has to evaluate: 



t e -x74T g -aT 



X(M) = / 7 = — dr, i€l,t>0. (91) 

V47TT 



2.3 



First we observe that: 



x(o,t) 



t e -ar 



-At 



1 2 



\fai 



e V dy 



2A 



lo \/Attt 2 A y/n J 

after the change of variables y = \far. Because Erf(— u) = — Erf(it) we can write: 

X {0,t) = 7^ { Erf(^) - Erf(-A*)} . 

Now, for any iel: 



Erf(Vat) 



Vle^ Erff- 
A I \2a 



because: 



4A 
_r/_ 



x 
2~Vt 



'of 



3inh(|a;|-A), 



(92) 



(93) 
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Moreover, because Erf(±oo) = ±1, we haveEq. [93]), which actually reduces to Eq. (|92|l when .t = 0. Therefore, 
the fundamental solution of Eq. f81j) is: 

/(*, t) = e- at G(x, t) + ^ |e«VS Erf f _£_ + ^ - e"^ Erf ( - y^tH - ^ sinh(M^), (94) 

which can be rewritten as Eq. (j88|) . □ 
Remark 7.5. With the choice (|86|) the process: 

B(i(t)) = S(M+(l-6t)i(t)), * > 0, 
has marginal density (|88| . Observe that, for any t > 0: 

B(M + (1 - 6t)j(t)) = b t B(t) + (1 - b t )B{j{t)), 
which naturally corresponds to Eq. [88]) - 

Remark 7.6. We observe that Eq. (j88| reduces to G(x,t) when a = (see Figure [8]), which is as expected 
because the memory kernel disappears. For small times, the nondocal memory effects are negligible and the 
process appears Markovian. Fig. 8 displays the fundamental solution at fixed t and various values of a, whereas 
Fig. 9 displays the fundamental solution at fixed a and various values of t. 



For large times we have: 
where: 



lim f(x,t) = lim <f>(x,t) = <j)(x), 

t — >oo t — >oo 



0(ar) = ^(cosh(aV^) - sinh(|a;|^)) - ^e"!^, x G 



(95) 
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Figure 9: Plot of the fundamental solution f(x,t), Eq. (|88|. at time t = [0.5, 1, 1.5], and a = 1. The dashed 
line represents the asymptotic distribution (j)(x), Eq. (|95|) . 

Remark 7.7. In view of Eq. ([82|) . it is always possible to choose the random time process l(t), t > 0, such 
that it becomes stationary at large times, in the sense of finite-dimensional densities. With this choice, the 
subordinated process B(l(t)) tends to a stationary process with asymptotic marginal distribution given by Eq. 
I|95p . For instance, if we look at Eq. lj87j) . as t — > oo we have l(t) = X ~ ae~ ar , t > 0. A less trivial example 
can be constructed by replacing the random variable X with a stationary process X(t), t > 0, such that for 
each t > the random variable X(t) has an exponential distribution with mean E(X(t j) = a -1 . The resulting 
process l(t) = min(X(£),i) is not increasing, has marginal distribution defined by Eq. (|82|) and tends to X{t) 
for large t. See Fig. 10. 

Remark 7.8. To obtain an idea on how fast the stationary regime is reached, one can look at the variance of 
the subordinated process. Using Eq. f83|) with m — 1, we find: 



which, for large times, tends exponentially to 2/a (i.e. the variance of eq. [95|. 
7.4 Exponential-decay kernel with logarithmic scaling time 

What happens if we choose an exponential kernel K(t) = e~ at and a logarithmic scaling time? That is: 



E{B{l{t)f) = -(l-e- at ), 



(96) 



g{t) = log(f + 1), t>0. 



(97) 



Since g'(t)K(g(t) - g(s)) = (t+ l)~ a (s + 1) 



, we get: 




(98) 
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a=l 



B(l(t)) 




Figure 10: Trajectory of the process B(l{t)) (top panel), with < t < 10, E(B(1)) = 1, l(t) = min(i,X(t)) 
where is an exponential White Noise with mean one. The corresponding trajectory of the random time 
l(t) process is presented in the middle panel. The estimated variance is computed on a sample of dimension 
N = 500. The smooth black line in the bottom panel corresponds to cr 2 (t) given by Eq. ([96)) and the stationary 
value is lim a (t) = 1. 

t — >oo 



Its fundamental solution is: 



/(*,*) 



(*+!)' 



-G(x,log(i+l)) 



— sinh(|a:|Vo) 



- i Erf 



2Vlog(t+l) 



-xVH Erf 



2 > /log(* + 1) 



- A /olog(f + 1) 



(99) 



Remark 7.9. As in Remark \7.7\ consider a random time process l(t), t > 0, with marginal distribution defined 
by Eq. (j82")) . that becomes stationary for large times. The subordinated process B(l(\og(t + 1))), t > 0, has 
marginal density function defined by f(x,t) of Eq. ([99| . Observe that in this case the random time process 
l(log(t + 1)) is no longer asymptotically stationary. This is because the translational time-invariance is broken 
by the logarithmic transformation. However, we can always consider a random time process I (t), t > 0, with 
the same marginal distribution of I (log(i + 1)), which becomes stationary for large times. Thus, the process 
B(l (t)) still has a marginal density function defined by f(x,t) but becomes stationary as t — > oo, in the sense 
of finite-dimensional distribution, with asymptotic marginal distribution given by Eq. (j95| . See Fig. 11. 



Remark 7.10. While B(l(t)), t > 0, satisfies Eq. f96|) and thus has a variance which tends exponentially fast to 
the limit value 2/ a, here the stationary regime is reached more slowly. Indeed, the variance of the subordinated 
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process is: 



E(B(lW 



- 1 



a \ (t + l) a J ' 

which, for large times, converges to the stationary value 2/a with a power-like behavior. 



(100) 
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Figure 11: Trajectory of the process B(l(log(t + 1))) (top panel), with < t < 10, E(B(1)) = 1, l(t) = 
min(i, X(t)) where X(t) is an exponential White Noise with mean one. The corresponding trajectory of the 
random time process Z(log(f + 1)) is presented in the middle panel. The estimated variance is computed on 
a sample of dimension N = 500. The smooth black line in the bottom panel corresponds to Eq. I|96p . The 
stationary value is lim <7 2 (log(£ + 1)) = 1. The stationary regime is achieved more slowly than in the case of 

t — >oo 

Figure E3 



8 Examples involving other diffusions 

We shall now consider examples of fractional and stretched Fokker-Planck equations involving diffusion operators 
other than V x = d xx which corresponds to standard Brownian motion. We will choose K{t) = t 13 ^ 1 /T((3) and 
g(t) = t a l® as in Section [Tl 1 and consider the partial integro-differential equation: 

u{x,t) = u {x) + —-- J av- 1 ^-* - 8 V) V x u(x lS )ds, < /3 < 1, a>0. (101) 

Its fundamental solution is the marginal density of the process: 

V(t) = Q{l {t a ^)), (102) 
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where Q(t), t > 0, is the stochastic diffusion associated to V x and lp(t), t > 0, is a suitable self-similar random 
time process. One has the following particular cases: 

• When a = and < < 1, Eq. (|101f) becomes the "time-fractional" Fokker-Planck equation: 

1 /"* 

u(x,t)=u (x) + w ^r (* - sf' 1 V x u(x, s)ds, < < 1, (103) 

r(p) 7o 

whose fundamental solutions are the marginal distributions of the process: 

V{t) = Q(lp(t)), 

and are given by: 

fv(x,t)= / f Q (x,r)f h (T,t)dT, (104) 



where 

f h (T,t)=t-' 3 M {Tt-P), T,t>0, O<0<1 (105) 
and /q(x, i) is the probability density of Q(t). 

• When = 1 and a > we get a "time- stretched" Fokker-Planck equation: 

u(x,t) = u (x) + [ as a - 1 V x u{x,s)ds, < < 1. (106) 
Jo 

In this case /;(r, i) = <5(t — t Q ) and we get: 

Mx,t) = f Q (x,t a ), 
which corresponds to the "stretched" diffusion: 

V(t) = Q(t a ), a > 0. 

• The case a = = 1 is trivial and corresponds merely to the Markovian case where the equation is: 

u(x, t) = uq(x) + / V x u(x,s)ds 
Jo 

whose fundamental solution is the density function of D(t) = Q{t), namely the Markovian process. 

In the following subsections we study the above equations under particular choices of the Fokker-Planck 
operator V x . In all the cases considered, we also give the results involving the exponential-decay kernel. 

8.1 Brownian motion with drift 

Let fj, £ R and a 2 > be given. Consider a linear diffusion = £?(^< <T ) on R satisfying the stochastic 

differential equation: 

dB^\t) = ndt + adB(t), t > 0, (107) 

where B(t), t > 0, is a "standard" Brownian motion. The process B^'(t), t > 0, is called Brownian motion with 
drift [i. It corresponds merely to a Brownian motion plus a drift term, namely: 

B M (t) = fxt + aB(t), t>0. (108) 

The marginal density function of B^>(t), t > 0, is: 

f BM (x, t)= . . 1 — exp (- (:E ~ f ) ) , t>0, (109) 



|CT|V4^rf V ° 4t 
which is the fundamental solution of the Fokker-Planck equation: 

d t u(x, t) = — fid x u(x, t) + <r 2 d xx u(x, t), t>0. (110) 
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8.1.1 The /3-power kernel. 

We consider the "fractional" Fokker- Planck equation, see Eq. (|103p (see also [35]): 

1 f* 

u(x,t) = u (x) + —- {t- s)^ 1 {-^d x u{x,s) + o 2 d xx u{x,s))ds, < p < 1. (Ill) 
1 IP; Jo 

Its fundamental solution can be regarded as the marginal density function of the process: 

D(t) = BV>(l fi (t)), t>0, Q<0<1, (112) 

where the process lp{t), t > 0, is a self-similar random time process with parameter H = [3/2, independent of 
B^', such that its marginal distribution is given by Eq. I|105p . 



Proposition 8.1. The fundamental solution of Eq. hill}) is: 

poo 

f D (x,t) = f BM (x,T)fl j3 (T,t)dT, t>0, XER, 



1 / (x~^t) 2 



f D {x,t)= / cxp - v ^ ' \M {T,t)dT, t>0, xel, (113) 

Jo HV47TT V 4(7 T / 

which is equal to: 

fn(xt ) 1 f- t^ll^ t -m H ^(\ X(J -M t -m (l/2,l/2),(l-/3/2 + /3fc,/3/2) \ 

J D {x,t)-e 2\a\^ k\ a 2,*\\x° \ l (0, 1), (fc + 1/2, 1/2) )' 1 j 

where the Fox H -function is defined by Eq. $67\). 
Proof: In order to evaluate fo(x,t) we write: 

/•OO 

f D (x, t) = lo-rh^'/ 2 " / e-^ T ' ia 'G(x', t)M p (t, t)dr, 
Jo 

where G(x,t) is the standard Gaussian density, see Eq. |(9]) and x' = xja. In view of Eq. ijTTj) . we have to 
evaluate an integral of the form: 

1 f°° 

®(x,t) = - e- aT M 1 / 2 (\x\,T)Mp{T,t)dT, xeR, t > 0, a > 0. (115) 

*(x,t) = - / e- a -T- l l 2 M l/2 {\x\T- l l 2 )t-PMp{Tt-e)dT 
* Jo 

= 5 jf \ Ml ' 2 (j) iye- ay2 t- M (y 2 t- )dy- 

after the change of variables y = \pr. Because of the symmetry, it is enough to consider only the case x > 0. 
We get: 



One has: 



where 



$(x,t) = ^(M 1/2 *Y t )(x), x>0, 



(<p*</>)(x)= I -<p ( -) (j)(y)dy 

Jo y \yJ 



indicates the Mellin convolution and where: 

Y t (x) = 2xe- ax2 t- M f3 {x 2 t-P), x>0, t > 0. (116) 
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Figure 12: Plot of the fundamental solution f(x, t), Eq. (|114p . with j3 = 1/2, at time t=l, for different values 
of the parameters fx = [1, 1.5, 2] and cr = [1, 1.5, 2]. 



Using the Mellin convolution theorem we get: 

M{2$(x, t);x,u] = M{M 1/2 {x);x,u}M{Y t (x);x,u}. (117) 
Because ofEq. [68J) andEq. [67J), this can be written as: 

M{2$(x, t);x, u} = r(1/ ^y , 2) M{Y t (x);x, u). (118) 

We now evaluate: 

fOO 

M{Y t (x);x,u} = / e~ ax2 2xt-l 3 M l3 {x 2 t-P)x u - 1 dx. 







After the change of variables x 2 t 13 = z, we get 

/>oo 

M{y t (a;);a;,ii}= / {ztP)^ u -^ e~ aztP M p {z)dz 
Jo 



= t #(„-i)^(z|J_ [°° z *-i+f M p(z)dz 
k=o ' ^° 

** Z! T, M{Mg(a:);a:, fc + 1/2 + u/2} 



!(„_!) f> (-a^) fc T(l/2 + k + u/2) 

^ fc! r(l + Pk - (3/2 + (3u/2) ' 
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Figure 13: Plot of the fundamental solution f(x, t), Eq. I|114p with /3 = 1/2, = 1, a = 1, at times t = [0.1, 1,2]. 



where we have used Eq. (|68|) and Eq. (j67|) . Thus: 



M{2$(x,t);x,w} = ^ 



(-at")* 



r(ti)r(i/2 + fc + tt/2) 



k\ T(l/2 + u/2)T(l + /3k- P/2 + f3u/2) 

Inverting the Mellin transform, Eq. [67| gives: 



k=0 



k=0 



(l/2,l/2),(l-/3/2 + /3fc,/?/2) 
(0,1), (A; + 1/2, 1/2) 



with x £ R and t > 0. Therefore, the fundamental solution of Eq. I|103p can be expressed as: 



)ix/2a 2 



^71 51 



-ftP/A^ 1 )* 



2|cr| ^ fe! 

1 1 fc=0 



(1/2, 1/2), (1-/3/2 + /%, 0/2) 
(0,1), (A +1/2, 1/2) 



that is Eq. ^THJ). □ 

When = and cr = l,Eq. 1114ft reduces to: 

f n ( x t ) Ir^H 2 ' (\x\t-W (1/2 ' 1/2) ' (1 ~ 0/2 ' m ) 
f D[x ,t)- 2 t H 2 2 y\x\t (0,1), (1/2, 1/2) )' 

that is, using the reduction formula for the Fox .//-function [30], 



f D {x,t) = \t-^Hl$(\x\t-W 



(l-/3/2,/3/2) 
(0,1) 



-Mp /2 (\x\,t). 



(119) 



31 




Figure 14: Trajectory of the process D(t) = B^> (h/ 2 (t)) defined in Eq. ITT2"|) with (3 = 1/2 (top panel). The 
random time process is l\/-z(t) = \b(t)\, where b(t) is a "standard" Brownian motion (middle panel). The variance 
and the mean are evaluated over a sample of size N = 5 ■ 10 4 and fit the theoretical values (bottom panel). 



As expected, we recover in this case the fundamental solution of the time-fractional diffusion equation (|72|) . 
Moreover, if we set j3 = 1 in Eq. (| 1 14|) we have: 



f D (x,t) = e^ 2 ° 



2 a 



k=0 



(1/2, 1/2), (1/2 +k, 1/2) 
(0,1), (1/2 + A, 1/2) 



E 

fc=0 



{-ixHj^f 1 , 



2\a 



(1/2,1/2) 
(0,1) 



1 



|cr|V47ri 



exp 



4a 2 t 



and we recover fsM (x,t). 



In Figure mi and Figure 1X51 we have used Eq. I|113p to plot the fundamental solution (|1 14f) with [3 = 1/2 
for different values of the parameters fi and a at fixed time and, for fixed parameters, at different times t. As 
expected, the fundamental solution is not symmetric in space with a time-growing skewness. Moreover, due 
to the presence of the positively taken drift term (// = 1), the probability to find the particle in the positive 
semi-axis increases with time (fig. [T3|) . 



In Figure fl4l is presented a trajectory of the process D(t) = B^{lp{t)) with (3 = 1/2. Using Eq. fT3|) it is 
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Figure 15: Marginal density function f(x,t) of the process B M (k/ 2 (t)) at time t = lEq. HTty . The histo gram 
is evaluated over TV = 10 5 simulated trajectories (see Figure fM)) . 



easy to write all the moments of the process: 



E(D(ty 



[m/2] 

E 

fc=0 



m\ 2k\{m- k)\ 2k m _ 2k 
■a \i 



2k 



t f3(m-k) 



T((3{m - k) + 1) ' 



</3 < 1, 



where m is an integer greater than zero and [a] indicates the integer part of a. Therefore, we have: 



m{t) = E(D{t)) = n 



and 



a 2 {t) = E{D{tf ) - m{tf = 2/i 2 



2a 2 



(120) 



(121) 



(122) 



r(2/? + i) r(/3 + i) 2 1 " r(/9 + i)' 

In the bottom panel of Figure E] the mean and the variance have been estimated from a sample of trajectories 
of the process B^\li/ 2 ){t). Then, they have been compared with the theoretical values given above. In Figure 
HU we compare the theoretical density function f(x,t) given by Eq. I114[) at time t = 1 with an histogram 
evaluated over a sample of N = 10 5 trajectories. 

8.1.2 Exponential-decay kernel 

The exponential- decay kernel case is straightforward. The non-Markovian Fokker-Planck equation is: 

ri 



i(x, t) = uq(x) + / e s \-fid x u(x, s) + a 2 d xx u(x, s))ds, a > 0. 
Jo 



(123) 
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If we indicate by Q(x,t) the fundamental solution of the Markovian equation; i.e. Eq. (|109|) 

Q{x, t)= ) cxp (- ^ ~ f) V t > 0, x £ R, (124) 



|CT|V47rf V 4fj2t 

then, using Eq. i|82|) . the fundamental solution of Eq. (( 1 23[) is: 

f(x, t) = e~ at g{x, t) + (1 - e- at )$(:c, f), (125) 

where: 

$( x t ) = - er x/2a / - - dr 

' l-e-«* Jo |a|V4^ 

Using Eq. ([H]) and Eq. ([93]) we have: 

Proposition 8.2. The fundamental solution of Eq. U23]) is: 



a6 !^l_ sinh (\xa- 1 \Ja+^ 



j -hV' + £ Erf fe + V(-^ 



cxp -rV & + 4& Erf yb-v( a+ ^> • (126) 



When i — > oo we obtain the stationary distribution: 



'2a 2 - \xa- 1 \\/a+ £^ j . (127) 



that is: 

<?0) = 7= 

2kl\/ a +4S 

8.2 Geometric Brownian motion 

Let /i e K and a 2 > be given. Consider a linear diffusion SonR defined by the stochastic differential equation: 

dS(t) = fj,S(t)dt + aS{t)dB{t), t > 0, (128) 

where B(t), t > 0, is a "standard" Brownian motion. The process S(t), t > 0, is called Geometric Brownian 
motion. If 5 starts in xq at time t = (i.e. P(S(0) = Xo) = 1), then a solution of Eq. I|128p is: 

S(t) = x exp [(n - a 2 /2)t + aB(t)] , t>0, x > 0. (129) 

The marginal density function of S(t) is the log-normal distribution: 

2 > 

1 



log(x/x )- (fi-a 2 /2)t\ 
fs(x,t)= ; exp I 5- ^] jt > ,x>0. (130) 

The function fs(x,t) is a solution of the Fokker-Planck equation: 

d t u(x, t) = [{2a 2 - fi) + (4a 2 - fi)xd x + a 2 x 2 d xx ] u{x, t), x > 0, (131) 
with deterministic initial condition 

uq (x) = S(x — 2:0)7 x>0, xq > 0. (132) 
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8.2.1 /3-power kernel 

If we introduce the /3-power kernel K(t) = r^)" 1 ^" 1 , < < 1, in this setting we obtain the following 
"fractional" Fokker-Planck equation: 

t) = u (a;) + — — (t- sf' 1 [{2a 2 - fj,) + (4a 2 - y)xd x + a 2 x 2 d xx ] u(x, s)ds, x>0. (133) 
1 IPj Jo 

A solution of the above equation with initial condition given by Eq. (|132fl is given by (see Corollary 16. ip : 



f D ( X ,t)= f S (x,T)Mf3(T,t)dT (134) 

Jo 

which is the marginal distribution of the process 

D(t) = S(l (t)), t>0, 0</3<l, (135) 

starting almost surely in xq > 0, where lp{t), t > 0, is a self-similar random time process with H = (3/2, 
independent of the geometric Brownian motion S(t) and with marginal density function given by Eq. (|105p . It 
is easy to see that: 

lD\x,t) = exp I ^ J-y e Mi/ 2 (F | , t)JW^ (t, i) dr, 

where: 

a=(ji- a 2 /2) 2 /4a 2 , x' = \og{x/x Q )/a. 
We have the same integral as in Eq. (|115p . Therefore: 
Proposition 8.3. for each t > 0; 

x\a\ \ 4cH 

^ i / (/i _ (T 2 /2)2 ^ y / (l/2,l/2),(l-/3/2 + /3fc,/?/2) \ 

x Z.fc?l v 4cr 2 J * «2,3^l* (0,1), (A + 1/2, 1/2) J' ( " fa) 

This result can be obtained directly from Eq. (|114p because our process is: 

D(t)=x eM Bitl 'Hb(t))), 

where > is a Brownian motion with drift p! = (fi — c 2 /2). When (3 = 1 we recover Eq. (j 130[) . Moreover, if 
fx = a 2 /2 (i.e. // = 0) we have (see Figure fT6|) : 



f D {x,t) = -^-t-^ 2 M p/2 



\og{x/x Q ) 



t~ p/2 ) , x>0, t>0, (137) 



x\a 

which is the marginal probability density of: 

D(t) = x e aBil ^ t) \ t > 0. 

In Figure [16] we show the plot of the fundamental solution f(x,t) in the particular case given by Eq. (|137p . 
Here we can see the behavior of the solution varying the parameter (3. For (3 = 1 we recover the log-normal 
density (|130|) with fj, = a 2 /2. In Figure [T7| we point out the dependence of the solution with respect to the drift 
parameter /i for fixed /3 = 1/2, t = 1, a = 1 and x = 1. In Figure fT8l we present the time evolution of the 
fundamental solution with (3=1/2 and (3 = 1/4. 
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Figure 16: Plot of the fundamental solution f(x,t), Eq. (|136j) at time t = 1, when /i = <r 2 /2 = lEq. I137p . 
x = 1, for different values of the parameter (3 = [1/4, 1/2, 3/4, 1]. For /3 = 1 f(x, t) reduces to the log-normal 
density (| 130[) . The angular point corresponds to the initial value xq = 1 and is due to the presence of | log(cc/a;o)| 
in the solution. 

In Figure[l9]we present a trajectory of the process D(t) = S(lp(t)) with (3 = 1/2, Eq. I|135p . We shall now 
compute the mean and the variance of the process D(t). We have that: 



E(S(t)) = E (x exp [(n - a 2 /2)t + oB{t))] ) = x Q cxp [(> - cr 2 /2)t] E 
Therefore, because: 

E(e" B V)= ( e° x G(x,t)dx = e° 2t -^= ( e' 1 ^ 1 dx, 
Jk v 47rt J« 



we have: 

E(S(t)) = x exp [(/x + a 2 /2)t] . (138) 

In the same way one has: 

E{S{t) 2 ) = x 2 cxp [{2^1 + 3a 2 )t] . (139) 

Using the above equations we have: 

E(S(l (t))) = x E (exp [(p + a 2 /2)l p (t)] ) = x £ ^ + ° 1 ' E(lp(t) k ), 



k=0 



which, using Eq. i(73)l . becomes: 



E(S(l (t)j) = x jr = ^((/i + ^ 2 /2)^), 



k=0 
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Figure 17: Plot of the fundamental solution f(x, t), Eq. (|136fl . with f3 = 1/2, a = 1, Xo = 1, at time t = 1, for 
different values of the parameter /U = [0.1, 1/2, 1, 2]. For fj, = 1/2 we have Eq. (|137[) . see also Figure fl6l 



where Ep(z) = V] z k /T((3k + 1) is the Mittag-LefHer function of order (i [26]. Similarly: 

fc=0 

^S^)) 2 ) = z 2 ^ (exp [(20 + 3a 2 )^(t)] ) = ^^((2/i + 3<r V)- 

Finally one has: 

m(t) = E(D(t)) = x E p ((p + a 2 /2)tP) 

a 2 {t) = E(D(t) 2 ) - m(t) 2 = x 2 [E fJ ((2^ + 3a 2 )t> 3 ) - Ep{& + a 2 /2)t?) 2 ] 
8.2.2 Exponential-decay kernel 

We now consider the exponential- decay kernel K(t) = e~ at , a > 0. The non-Markovian Fokker-Planck equation 
is: t 

i(x,t)=u (x) + e~ a{t - s) [(2a 2 - Li) + {4a 2 - fi)xd x +<7 2 x 2 d xx ]u{x,s)ds, a > 0. (141) 



(140) 



o 



We denote by Q(x, t) the fundamental solution of the Markovian equation; namely Eq. (|130[) 

/ , 



G&'t) = i //I— 7 ex P 



V 



log(x/a;o)-(//-cr 2 /2)i)' 
^44 



, ic,t>0. (142) 



Then, using Eq. (j82|) . the fundamental solution of Eq. (|123fl is: 

f(x, t) = e- at G(x, t) + (l- e~ at )<S>(x, t), (143) 
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where: 



and: 



$0M) = T 



a 
\a\x 



e log(^)(^-) 



G(x',T)e- aT dr 



, x > 



, _ (ti-a 2 /2) 2 



Thus, as in Eq. (|94|) . we have: 
Proposition 8.4. 



4(T 2 i 4a, a,/i>0, cr > 0, 
x' = log(x/a;o)/(7, a; > 0. 



f( X ,t)=e- at G(x,t) 



■ ^i^C-^S 1 ) J _L_ ex p(a;'V^)Erf ( + V^t ) - cxp(-xVa')Erf 



4^ 



4|cr|xVa' 

The stationary distribution, obtained as t — > oo, is 



e log( ^ )(li ^ ) sinh(|a;'|Va 7 ). 



= „i i° r~7 CX P ( lo S f — ) ( 



a; \ / 2/i — (j^ 
, V V 4a 2 



(cosh(a;'v a') — sinh(|a/|va') 



(144) 



(145) 



(146) 



38 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 19: Trajectory of the process D{t) = S(ly 2 {t)) defined in Eq. ljT53)l with /3 = 1/2 (top panel). The 
random time process is £1/2 (t) = \b(t)\ (middle panel). The variance and the mean are evaluated over a sample 
of size N = 5 • 10 4 and are presented together with the theoretical functions, Eq. (| 140|) . in the bottom panels. 



9 Conclusions 

Theorem 16 . II states that the fundamental solution f{x, t) of a non-Markovian diffusion equation of the form (f43| 

u(x,t) = u (t) + [ g'(s)K (g(t) - g(s))V x u(x, s)ds, t > 0, (147) 

Jo 



f(x,t) 



G( X ,T)h(T,g(t))dT, 



(148) 



where Q{x,t) is the fundamental solution of the Markovian equation l(4T|) and h(r, t) is the fundamental solution 
of the non-Markovian forward drift equation 

u(r,t) = u (t) - ( K(t- s)d T u{T,s)ds, T,t>0, (149) 
Jo 

If the memory kernel K(t) is chosen in a suitable way (see Section [3]), the solution f(-,t) preserves non- 
negativity and normalization for all t > 0. Thus, it can be interpreted as the marginal density function of a 
non-Markovian stochastic process. In view of Eq. (| 148|1 . this stochastic process is naturally interpreted as a 
subordinated process Eq. (|33)l . 



We focused on two kind of memory kernels: the power kernel K{t) = t 13 " 1 ^^), < /3 < 1, and the 
exponential-decay kernel K (t) = e~ at , a > 0. 
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The first provides the so-called time-fractional Fokker-Planck equations (|101[) . In particular we studied 
the case V x = d xx (see Section 17. ip , which corresponds to the choice of a "standard" Brownian motion for the 
parent Markov model. In this case, the fundamental solution can be written in terms of an entire transcendental 
function, see Eq. iff^l . and is related to a Fox i/-function through Eq. ([68| . We have also considered more 
complicated cases, namely Brownian motion with drift fi (see Section 18. ip and Geometric Brownian motion 
(see Section I8.2|) . In these cases the fundamental solutions can be written in terms of a superposition of Fox 
^-functions, see Eq. fTTiland Eq. ((Tgg)! . 

The exponential-decay kernel corresponds heuristically to a system in which the non-local memory effects 
are negligible for small times. In fact, the fundamental solution can always be written in the form of Eq. (|88l) . 

f(x, t) = e- at G(x, t) + (1 - e- at )(j)(x, t), t > 0, 

where Q(x, t) is the fundamental solution of the Markovian equation, and where the function (j)(x, t) is a proba- 
bility density which becomes stationary for large times. Therefore, it is always possible to find stochastic models 
that become stationary for large times and whose marginal density is given by Eq. (| 1-48|) . 

However, see Subsection the stochastic representation is not unique, that is, there are many different 
stochastic processes whose marginal density is f(x,t). For example, consider the case where V x = d xx and 
g(t) = t. Then, /(x, t) is the marginal density of B(l(t)), t > 0, where B(t) is a"standard" Brownian motion and 
where l(t) is a random time process satisfying Eq. (| 149|) . If the random time l(t) is required to be self-similar 
of order /3, then in view of Theorem l3.il the memory kernel must be a power function K(t) = t ~ l /T(/3) with 
< P < 1. The corresponding non-Markovian diffusion equation l|147p is called in this case time-fractional diffu- 
sion equation of order j3 (see Subsection l7.ip . The corresponding random time process l(t) = lp{t), can be the lo- 
cal time of a d = 2(1 — /3)-dimensional fractional Bessel process or, alternatively, the inverse of the totally skewed 
strictly /3-stable process. However, f{x,t) is also the marginal density of the process Y(t) = y/lp{\)Bp/ 2 (t), 
where Bp/ 2 ls a fractional Brownian motion independent of the random time lp(t). In all the previous exam- 
ples, the self-similarity parameter H = j3/2 is restricted to the region < H < 1/2. We can obtain stochastic 
processes with higher values of the self-similarity parameter by introducing the time-scaling function g(t). In 
this way, for example choosing g(t) = W' 3 , < a < 2, we obtain the process D(t) = _B(Z ) 3(£ Q / /3 )), t > 0, and 
the process y(t) = y/lp(l)B a / 2 (t), t > 0, which are self-similar with parameter H = a/2 so that < H < 1 
(see Subsection I7.2p . In contrast to D(t) the process y(t) has stationary increments. 

The solution of the "non-Markovian" equation (|147p can be stated explicitly in all the cases considered. We 
computed it analytically and graphed it in particular cases. This solution is a marginal (one-point) density 
function. We have then presented various random processes whose marginal density function coincides with 
that solution. 
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